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Abstract 

We present new third- and fifth-order Godunov-type central schemes for ap- 
proximating solutions of the Hamilton- Jacobi (HJ) equation in an arbitrary num- 
ber of space dimensions These are the first central schemes for approximating 
solutions of the HJ equations with an order of accuracy that is greater than two 
In two space dimensions we present two versions for the third-order scheme- one 
scheme that is based on a genuinely two-dimensional Central WENO reconstruc- 
tion, and another scheme that is based on a simpler dimension- by-dimension re- 
construction The simpler dimension-by-dimension variant is then extended to a 
multi-dimensional fifth-order scheme Our numerical examples m one, two and 
three space dimensions verify the expected order of accuracy of the schemes 
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1 Introduction 

We are interested in high-order numerical approximations for the solution of multi- 
dimensional Hamilton-Jacobi (HJ) equations of the form 

4>t + H(\7<f>) = 0, x=(xi, Xd)eR d , 

where H is the Hamiltonian, which we assume depends on V</> and possibly on x and t 
In recent years, the H J equations attracted a lot of attention from analysts and numerical 
analysts due to the important role that they play in applications such as optimal control 
theory, image processing, geometric optics, differential games, calculus of variations, etc. 
The main difficulty m treating these equations is due to the discontinuous derivatives 
that develop in finite time even when the initial data is smooth Vanishing viscosity 
solutions provide a good tool for defining weak solutions when the Hamiltonian is convex 
- [15] The celebrated „ viscosity . solution provides a suitable extension of weak solutions 
for more general Hamiltonians [3, 7, 8, 9, 10, 28, 29]. 

Given the importance of the HJ equations, there has been relatively little activity 
developing numerical tools for approximating their solutions This is surprising given 
that most of the numerical ideas are based in the similarity between hyperbolic conser- 
vation laws and the HJ equations, and the field of numerical methods for conservation 
laws has been flourishing in recent years 

Converging first-order approximations were introduced by Souganidis in [38] High- 
order upwind methods were introduced by Osher, Sethian and Shu m [34, 35] These 
methods are based on Harten’s Essentially Non-Oscillatory (ENO) reconstruction [13], 
that is evolved in time with a first-order monotone flux The Weighted ENO (WENO) 
mterpolant of [18, 32] was used for constructing high-order upwind methods for the HJ 
equations m [17], and extensions of these methods for triangular meshes were introduced 
in [1, 40] We note in passing that there are other approaches for approximating solutions 
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of HJ equation such as discontinuous Galerkin methods [14, 24] and relaxation schemes 

[ 20 ] 

A different class of Godunov-type schemes for hyperbolic conservation laws, the so- 
called “central schemes”, have recently been applied to the HJ equations The prototype 
for these schemes is the Lax-Friedrichs scheme [11]. A second-order staggered central 
scheme was developed for conservation laws by Nessyahu and Tadmor in [33] The 
mam advantage of central schemes is their simplicity Since they do not require any 
(approximate) Riemann solvers, they are particularly suitable for approximating multi- 
dimensional systems of conservation laws. Lin and Tadmor applied these ideas to the HJ 
equations in [31] There, first- and second-order staggered schemes versions of [2, 19, 33] 
were written in one and two space dimensions An L 1 convergence of order one for this 
scheme was proved in [30]. After the introduction of a semi-discrete central scheme 
for hyperbolic conservation laws in [23], a second-order semi-discrete scheme for HJ 
equations was introduced by the same authors in [22] While less dissipative, this scheme 
requires the estimation of the local speed of propagation at every grid point, a task that 
is computationally intensive m particular with problems of high dimensionality By 
considering more precise information about the local speed of propagation, an even less 
dissipative scheme was generated in [21] 

Recently we introduced m [5] new and efficient central schemes for multi-dimensional 
HJ equations These non-oscillatory, non-staggered schemes were first- and second-order 
accurate and were designed to scale well with an increasing dimension Efficiency was 
obtained by carefully choosing the location of the evolution points and by using a one- 
dimensional projection step Avoiding staggering by adding an additional projection 
step is an idea which we already utilized in the framework of conservation laws [16] 

In this work we introduce third- and fifth-order accurate schemes for approximating 
solutions of multi-dimensional H J equations. These are the first central schemes for such 
equations of order greater than two This work is the HJ analog to the corresponding 
works in conservation laws an ENO based central scheme [4], and the Central WENO 
(CWENO) central schemes [25, 26, 27] We announced a preliminary version of the one 
dimensional results in a recent proceedings publication [6] 

The structure of this paper is as follows We start in §2 with the derivation of our 
one-dimensional schemes. A third-order WENO reconstruction scheme is presented in 
§2 2 This scheme required a fourth-order reconstruction of the point- values and a third- 
order reconstruction of the derivatives at the evolution points Even though the optimal 
location of the evolution points in one dimension is m the center of the interval, m order 
to prepare the grounds for the multi-dimensional schemes we write a reconstruction for 
an arbitrary location of the evolution points A fifth-order method is then presented m 
§2 3 

We turn to the multi-dimensional framework in §3 Here there is flexibility in the 
reconstruction step For simplicity we carry most of the discussion in two space di- 
mensions Extensions to more than two space dimensions are presented in §3 4 First, 
we provide a brief outline of the general structure of two-dimensional central schemes 
m §3 1 The mam remaining ingredient, the reconstruction step, is then described in 
the following two sections For a two-dimensional third-order scheme we present in 
§3 2 two ways to obtain a high-order reconstruction of the approximate solution at the 
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evolution points The first option m §3 2 1 is based on a genuinely two-dimensional re- 
construction An alternative dimension-by-dimension approach is based on a sequence 
of one- dimensional reconstructions and is presented m §3 2 2 Our numerical results 
show that both approaches are essentially equivalent Hence, the rest of the paper deals 
with the dimension-by-dimension reconstruction A fifth-order dimension-by-dimension 
extension of the one-dimensional scheme in §2.3 to two dimensions is then presented in 
§3 3 Since the solution at the next time step is computed at grid points that are differ- 
ent from those on which the data is given, we reproject the evolved solution back onto 
the original grid points Different ways to approach this reprojection step are discussed 
in §3 2 3. 

We conclude in §4 with several numerical examples in one, two and three space 
dimensions that confirm the expected order of accuracy and the high- resolution nature 
of our scheme We compare our results with the scheme of Jiang and Peng in [17] We 
also study the convergence rate after the emergence of the discontinuities m the solution. 

Acknowledgment: We would like to thank Volker Elling for helpful discussions through- 
out the early stages of this work The work of D. Levy was supported in part by the 
National Science Foundation under Career Grant No DMS-0133511 


2 One-Dimensional Schemes 

2.1 One-Dimensional Central Schemes 

Consider the one-dimensional Hamilton-Jacobi equation of the form 

4>t{x, t) + H (4> x ) = 0, x € R (2 1) 

We are interested in approximating solutions of (2 1) subject to the initial data <j>(x,t = 
0) = 0o (%) For simplicity we assume a uniform grid grid in space and time with mesh 
spacmgs, Ax and At, respectively Denote the grid points by x t = lAx, t n = nAt, and 
the fixed mesh ratio by A =-At/Ax Let denote the approximate value of_0 (x t , t n ), 
and (<£*)" denote the approximate value of the derivative <p x {x u t n ). We define the 
forward and backward differencing as A + <p" = tp* +1 — </>" and A = <p* — 

Assume that the approximate solution at time t n , is given A Godunov- type 
scheme for approximating the solution of (2 1) starts with a continuous piecewise- 
polynomial <p(x,t n ) that is reconstructed from the data, <p?, 

(p{x, t n ) = ^ P l+ 1 (x, t n )x l+ i 2 (x) (2 2) 

t 

Here, ^+ 1 / 2 ( 2 ) is the characteristic function of the interval [x*,x 1+ i], and P»+i/2(^i t n ) 
is a polynomial of a suitable degree that satisfies the interpolation requirements 


±( x t+{3)t ) — — 0, 1 
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The reconstruction (2 2) is then evolved from time t n to time t n+1 according to (2 1), 
and sampled at the half-integer grid-points, {x t+ i/ 2 }, where the reconstruction is smooth 
(as long as the CFL condition A \H' (<p x )\ < 1/2 is satisfied) 


«t n+1 

= <4 - l » {v, (*.+} , t )) ir (2 3 ) 

The point- value V^I+i/2 * s obtained by sampling (2.2), at x, +1 / 2 , i e. <^" +1 / 2 = <p(x t+ i/ 2 , t n ) 
Since the evolution step (2 3) is done at points where the solution is smooth, we can ap- 
proximate the time integral at the RHS of (2 3) using a sufficiently accurate quadrature 
rule For example, for a third- and fourth-order method, this integral can be replaced 
by a Simpson’s quadrature, 


pt n + 1 A i 

/, ^(^(^.-))*- T [»(« i )+«(<d)+^(cr)] < 24 > 


The derivative at time t n , y >' % +\/2 ls obtained by sampling the derivative of the recon- 
struction (2 2), i e , V^+ 1/2 = t n ) The intermediate values of the derivative m 

time, y/ + T/f, and ^" 1 + / 1 2 , which are required m the quadrature (2 4) , can be predicted 
using a Taylor expansion or with a Runge-Kutta (RK) method Alternatively, (2 1) can 
be treated as a semi-discrete equation by replacing the spatial derivatives with their 
numerical approximations and integrating in time via an RK method. 

The only remaining ingredient to specify is the reconstruction (2 2) Below we 
present two reconstructions The first is a fourth-order reconstruction of the point- 
values and the derivatives which leads to a third-order scheme, and the second is a 
sixth-order reconstruction that results in a fifth-order scheme 


Remarks 

1. In order to return to the original grid, we project <p"_|^ 2 back onto the integer 
grid points {a:,} to end up with <p" +1 . This projection is accomplished with the 
same reconstruction used to approximate <p" +1/ , 2 fr° m V 5 / 

2 In order to maximize the size of the time-step, the evolution points should be" 

taken as far as possible from the singularities in the reconstructed 
piecewise-polynomial In one dimension the appropriate evolution point is 
located at x t+ i/ 2 In d-dimensions with a uniform grid with spacing Ax, the 
optimal evolution points are located at x l+a = x, + aAx m each direction, where 
OL — \/ ^ d (see [5]) One of the multi-dimensional schemes we present in 

§3 is based on one-dimensional reconstructions Hence, m order to prepare the 
grounds for the multi-dimensional setup, we write the one-dimensional 
reconstruction m this section assuming that the evolution points are x t ± a The 
reader should keep in mmd that in one dimension, a = 1/2. 

3 We would like to point out that one does not need to fully reconstruct the 
polynomials P, +1 / 2 (x, t n ) The only values that the scheme requires are the 
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approximated point- values <£>" +1 / 2 = v?(x,+i/ 2 , t n ) and the approximated 
derivatives i p' l+ ij 2 = fi( x %+i/ 2 ) Hence, m the rest of the paper whenever we refer 
to reconstruction steps we directly treat the recovery of these two quantities. 


2.2 A Third-Order Scheme 

A third-order scheme is generated by combining a third-order accurate ODE solver in 
time for predicting the intermediate values of the derivatives in (2 4), with a sufficiently 
high-order reconstruction m space 

Given m order to invoke (2 3), we should compute two quantities in every time 
step, the point- values at the evolution points, ip t ± a , and the derivative <^ ±Q . In order to 
obtain a third-order scheme, the approximations of the point- values should be fourth- 
order accurate, and the approximation of the derivatives should be third-order accurate. 
In this scheme, the reconstruction of the point-values is done in locations that are 
staggered with respect to the location of the data The reconstruction of the derivatives, 
which is required in every step of the ODE solver, is done at the same points where 
the data is given Since we anyhow need two types of reconstructions and due to 
symmetry considerations, we derive a fourth-order approximation of the derivatives 
Obviously, this more accurate reconstruction of the derivatives does not increase the 
order of accuracy of the scheme but it does reduce the error 


1 The reconstruction of <p t ± a from <p t 

A fourth-order reconstruction of <p t+a can be obtained by considering a convex 
combination of two quadratic polynomials, each of which requires the evaluation 
of on a three-point stencil One quadratic polynomial <p_(x) is constructed 
on a stencil that is left-biased with respect to x l+ctl {x t _i,x t . x, + i}, while the 
other polynomial <p+(x) is constructed on a right-biased stencil, { x„x t+ i,x t+ 2 }, 
see Figure 2 1. We set 


V’-.l+Q — 

V+,*+oc = 


= ( Q 2°‘ ) V>*-1 + (1 - a 2 ) <A + ) ¥>t+i , (2 5) 


2 3a + a 2 _ , 2 y , ( ~ 0L + a 2 

<p t + (2a -a ) <p t+1 + I — ) <p l+2 


For smooth <p, a straightforward computation shows that <p± it + a = <p(x t + a ) + 
O (Ax 3 ), and 


^ (2 - a) <p-, t+a + i (1 + a) <f +>t+Q = <P (x l+a ) + O (Ax 4 ) 

Similarly, the reconstruction of <p t ~ a is obtained using the quadratic polynomials 
<p_(x) based on the left- biased stencil enclosing x t _ Q , {x,_ 2 , x t _i, x,}, and <p+ (x) 
based on the right biased stencil x,,x t+ i}, 




—a + a 2 


Pi -2 + (2a — a 2 ) i + 


2 - 3a + a 2 


fft, (2.6) 
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Figure 2 1 The two mterpolants used for the third-order reconstruction at the evolution 
point at x l+a 




f a + a 2 

V 2 


) <p x - 1 + (l - a 2 ) <p x + 


f —a + a 2 

V - 


<A+ 1- 


This time, = <p (x,_ a ) + (9 (Ax 3 ), and 


| (1+ a) + § ( 2 - a ) P-h'-a = v 3 (®.-a) + O (Ax 4 ) 

A fourth-order WENO estimate of yp,± Q is therefore given by the convex combina- 
tion 


<A±a = w l±Q ip. rt±a + W+ ±a ip +>l±a) (2 7) 

where the weights satisfy w~± a + w+ ±Q = 1, wf± a > 0, V? In smooth regions we 
would like to satisfy w~ +a = w+_ a ~ (2 — o) /3 and = w~_ a « (1 + a) /3 
to attain an O (Ax 4 ) error When the stencil supporting <-p x ± a contains a discon- 
tinuity, the weight of the more oscillatory polynomial should vanish Following 
[18, 32], these requirements are met by setting 




E« 


a ;±« = 


<4a 


(e + S x ± a ) 




(2 8 ) 


where k,l € {+, — } The constants are independent of the grid index i and are 
given by c~ +a = c t t Q = (2 - a) /3, c+ a = c~_ a = (1 + a) /3 We choose e as 10 -6 
to prevent the denominator m (2 8) from vanishing, and set p = 2 (see [18]). The 
smoothness measures Sf should be large when ip is nearly singular Following 
[18], we take S t ± a to be the sum of the L 2 -norms of the derivatives on the stencil 
supporting <p± If we approximate the first derivative at x, by A + <p t /Ax, the 
second derivative by A + A _ <p,± Q /( Ax) 2 , and define the smoothness measuie 
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then we have S l+a = 5, [-1,0], S? +a = S x [0,1], 5,_ Q = 5, [-2, -1] and S+ a = 
S t [—1)0] 

For future reference we label the reconstruction in this section with the procedural 
form 


ip x ± a = reconstruct_<^_lD_3 («, ±a, ip) , (2 10) 

where <p is the one-dimensional array (ipi, ,<Pn)- This notation will be used in 
the dimension-by-dimension reconstructions in §3 

2 The reconstruction of <p\± a from <p x ± a 

The values of <p we recovered in the previous step at the regularly spaced lo- 
cations {x,± a } can be used to recover the derivative <p' x ± a via a (non-central) 
WENO reconstruction To obtain a fourth-order WENO approximation of <p' x ± a , 
we write a convex combination of three quadratic interpolants on the sten- 
cil {x x — 2ia j 3-d— liaj Po,t+a the StenCll l±a ; ^i+lia} 1 &nd 

on the stencil (x t ± a , x.+iia, x x + 2 ±a] For smooth ip, 

P—,x±a = 2 Ax ^Px— 1±<* d~ 3( p x ±a) = P (^lia) ~b O (Ax ) , 

Po,x±a = 2/S.x (P t+ ^ a ~ P*-^**) = P (®*±o) ~ O (Ax ) , (2.11) 

P+, i±a = 2Al^ ^Px+a d* &P x +l±a Px+1 ±a) = P (*^t±a) d“ O (Ax ) 

A straightforward computation yields 

\p'-,x±a + \p'o,x±c + £<4,*±a = p' (*t±a) + O (Ax 4 ) . 

- - . The fourth-order WENO estimate of p',± a from <p x ± a >s therefore 

Px±a = W 7± a p'-,x±* + dLaPo.ti* + W t±*P'+,x±«i ( 2 12 ) 

where the weights w are of the form (2 8) with k , l € {+, 0, — }, c~ = c + — 1/6, c° = 
2/3, and the oscillatory indicators are S~ ±a = S x ± a [—2, —1], S°± a = S x ± a [—1)0], 
and 5 t ± a = S x ± a [0, 1] 

For future reference we label the above reconstruction of <p' t±Q with the procedural 
form 


p[± a = reconstruct_</?'_lD_3 (i, ±a, <p ± a ) , 


(2.13) 


where tp± a is the one-dimensional array (pi± Q , 


<PN±a)- 
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We would like to summarize the one-dimensional third-order algorithm in the follow- 
ing, where RK (<p?± a , <p' t ± a > At) is the third-order Runge-Kutta method which integrates 
(2 1) and is used to predict the intermediate values of the derivatives Each internal 
step of the RK method will require additional reconstructions of (p' z ± a from that step’s 

<Pz±a 

Algorithm 2.1 Assume that {p™} are given 

1 Reconstruct- 

<Pi ±a = reconstruct-^ -ID .3 (i, ±a, c p n ) 
ip' t ± a = reconstruct -ip' AD .3 (i, ±a, pA± a ) 

t- 

2 Integrate 

= RK( V ; ±a ,^,At/2) 

p[2t 2 = reconstruct jp'.lD -3 

= RK {<P?± a , <p£ a , At) 
p >[ = reconstruct-ip' .ID .3 (i, ±a, ^±i) 

v-si = + y [«■ Wij + {v ,:: 1 ) + h ter 1 )' 

3 Reproject 

p™ +1 — reconstruct-ip .ID -3 (i, ^±i) • 

Remark It is possible to replace the Simpson’s quadrature in the integration step with 
a single RK time-step, <p™±l — RK (<£>” ±a , p' t ± a , At) . Our simulations show that this 
choice reduces the complexity of the computation but also reduces its accuracy 



10 


S. Bryson and D Levy 



Figure 2 2: The three interpolants used for the fifth-order reconstruction tp t + a at the 
evolution point at x t+Q In this example, because of the large gradient between x l+ \ 
and x t + 2 , the interpolant <^_ will have the strongest contribution to the CWENO recon- 
struction at x l+a . 

2.3 A Fifth- Order Scheme 

In order to obtain a fifth-order scheme, we need a sixth-order approximation of the 
pomt-values of <p, a fifth-order approximation of the derivative <p', and a higher-order 
prediction of the intermediate derivatives which appear m the quadrature formula Due 
to arguments similar to those given in §2 2, we again derive a more accurate reconstruc- 
tion of the derivatives, which in this case is sixth-order 

We start with the reconstruction of ip l+a from ip. We write sixth-order interpolants 
as a convex combination of three cubic interpolants, each of which requires the evaluation 
of <p on a four-point stencil We use the polynomials <p~(x) defined on the left-biased 
stencil {x,_ 2 ,x l -i,x,,x t +i} ) <po(x) defined on the centered stencil {x t _i,x„x, + i,x, +2 } 
and <p+(x) defined on the right-biased stencil {x,, x t+ i, Xj+ 2 , x t+ 3 }, see Figure 2 2 For 
smooth ip 

= a\<Pi -2 + i + a 3 ip t + a 4 <p t+1 = <p (x t+a ) + O (Ax 4 ) , (2 14) 

= astyt-l + 06^4 + ^7^1+1 + “8^1+2 = (^i+a) + O (Ax 4 ) , 

<P+,i+ a = &9<Pi + ai 0 ^ t +l + Gll<A +2 + ^12^+3 = <P (Xt+a) + O (Ax 4 ) , 

where the constants are given by 


" Ol 

1 1 3 

= -a— -or, 

-a 2 = 

1 O 1 Q 

—a + —a - + —o; - , 

6 6 

2 2 

a 3 

= 1 + — a 2 - 

1 3 

-2*’ 

1 1 2 1 3 
a4 ~3 a+ 2 a + 6 a ’ 


1 1 2 

1 3 

. 1 2 1 3 

a 5 

= — a H — a — 
3 2 

6*’ 

a 6 = 1 - -a - a + -a 


, 1 2 1 

3 

1 ,1 3 

dj 

= a + -a — — q 

' J 

as = ~g a+ g a = _ai > 

a$ 

= 1 - ^a + a 2 
6 

1 3 

or 

6 

, aio = 3a — -a 2 + -a 3 

a u 

3 2 

= --a + 2 or — 

Zj 

1 3 
2*’ 

1 1 2 1 3 
ai2 = 3 Q '2 a + 6 a 

At x t ~ a we have 




,1— or 


0i2^*-3 + an^i -2 + aio^i -1 + a 9 <Pt = <P (Xi-a) + O (Ax 4 ) , (2 15) 
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9?o,t-Q = <W *-2 + G7<A-i + a 6 <A + a 5 Wi = <p (£,_<*) + O (Ax 4 ) , 
<p+, t - a = a 4 ip t _i + a 3 ip x + a 2 ipr+i + aip t+2 = <p (x,_ a ) + O (Ax 4 ) 

A straightforward computation yields 


C i ± aP - S±a + C° l±a <f 0 , t ±a + l± a = ( Xi ± a ) + O (Ax 6 ) , 


where 




C f±a 


— s* 

'“'i+a 


+ 1 2 1 3 

C% ~ a — 2Q a 4 a+ 10’ 
1 2 1 3 

'l0 a + 10 a+ 5’ 

1,3 1 


C ~~ Q ~ 20° + ^0 a+ 10 


A sixth-order reconstruction of (p t ± a 1S therefore given by 


(2.16) 


<Pi ±a = W«toV , -,t±a + ^ t ±o^0,i±a + W,±aP+,t±Of> (2 17) 

where the weights w k are given by (2 8) with k,l € {+,0,—}, and the constants c k 
are given by (2 16) The oscillatory indicators are given via (2 9) by S~± a = S l [—2, 0], 
5 t ° ±a = 5, [-1, 1] and S& tt = 5, [0, 2] 

A sixth-order approximation of from ip t ± a is written as a convex combination 
of four cubic mterpolants This reconstruction is similar to the third-order case, and is 
based on a non-central WENO reconstruction We skip the details and summarize the 
result. 


PLc. = ™l±M,r± a + U&a^Kfca + W La^±a + ^±a^i±a» ( 2 - 18 ) 

where 

^l.tia = ( 2 ^ t— 3^° 9(^ 1 _2 ±q + H^±q)) 

^2,i±c* = (^*~^ a l±a "b 3c^ t ± a + 2 (^,^.i±q,), 

V^3,i±a = + 6<p t +i±o, ( ^ , t+2±a)) 

^4,«±a = + 18v?t+1±a “ 9< ^+ 2± “ + 2< ^+ 3±a) 

Here the weights in* are given by (2 8) with C\ = c 4 = 1/20, C 2 = c 3 = 9/20, = 

5 t±Q [-3, -1], S? ±a = S l±a [-2, 0], S?± a = S l±a [-1, 1] and S 4 ±Q = S t±Q [0, 2] 


Notations 


(2 19) 


1 We label the reconstruction of the point-values (2 17) as 
(p t ± a — reconstruct-^ _1D_5 (i, ±a, tp ) , 
where <p is the one-dimensional array , Pn) 
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2 We label the reconstruction of p>[± Q (2 18) as 

<p[± a = reconstruct_</ 7 '_lD _5 (i, ±a, <p± a ) , (2 20) 

where p>± a is the one-dimensional array (<pi± a , , Vn±o) 

Remarks 


1 To conclude, the fifth-order method is given by Algorithm 2 1 , where the 
fourth-order reconstructions are replaced by the sixth-order reconstructions 
(2 19)-(2.20) As is, this scheme is only fourth-order in time. A higher order 
method in time can be easily obtained by replacing Simpson’s quadrature with a 
more accurate quadrature and computing the sixth-order approximations for the 
point-values and the derivatives at the new quadrature points 


2 . We choose to predict the intermediate values of the derivatives m time using the 
fourth-order strong stability preserving (SSP) Runge-Kutta scheme of [ 12 ]. For 
s G { g, 1 } , the SSP-RK scheme is given by 


= 

<p( 2 ) = 

tpW = 


<P 


<p n --sAtH(<p n x ), 


5000 


+ 


n-fs 


649 n 10890423 n . , 951 ft) 

w H sAtH (<p2) H V' 1 ~ 

160(T 25193600 K ^ xJ 160(r 

53989 _ 102261 AirT . 4806213 m 

p> H sAtH ( (fiZ) H cp’ 1 ' 

250000CT 5000000 Wx} 20000000* 

7873 


7873 ^ W), 


5121 

20000 


+ 


sAtH (<p^) H <p^ -I sAtH (p>^) 
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“Alternatively, the Natural Continuous Extension of the RK method [39] can be _ 
used to produce the intermediate values p >' n+ 5 and p>' n+1 with a single RK step, 
though we observe that errors are somewhat larger in this case 


3 Multi-Dimensional Schemes 

3.1 Two-Dimensional Central Schemes 

Consider the two-dimensional HJ equation of the form 

4>t + H (V 0) = 0, x = (aq , X 2 ) G R 2 , (3 1) 

subject to the initial data cf>(x,t = 0) = <f>o{x) Denote x %t3 = (xi +iAxi,xi + jAx^) 
Similarly to the one-dimensional setup, <p tiJ will denote the approximation of at x ZJ 
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We define the two sets of grid points, 1+ = {x tiJ> x,+ij, and /_ = {x hJ ,x^ } , x hJ -i 

and denote by T+, T_ the triangles with vertices I + and /_ respectively For simplicity 
we assume a uniform grid Aar = Ax 2 — Ax 

Assume that the approximate solution at time t n , is given Similarly to the one- 
dimensional setup m §2.1, a Godunov- type scheme for approximating the solution of 
(3.1) starts with a continuous piecewise-polynomial (p{x,t n ) that is reconstructed from 
the data, , 

0(3 t n ) = Pjf (£, t n )xr ± {£) (3 2) 

As usual, Xt± (x) is the characteristic function of the triangle T±, and t n ) is a 

polynomial of a suitable degree that satisfies the interpolation requirements 

P^(x h t n ) = <p(xi,t n ), xi € /± 

(see Figure 3 1) The reconstruction (3 2) is then evolved from time t n to time t n+1 
by (3 1), and sampled at the evolution points {xi± atJ ±a} In two dimensions the choice 
a = 1/(2 + y/2) guarantees that the solution remains smooth at the evolution point as 
long as the CFL condition ^ | H' (V<p)| < a is satisfied The evolved solution now reads 

rt n+1 

tf&ia = V??±aj±a ~ / H ( V ^ (*.±oj±o, t)) A\ (3 3) 

The point-values <p" ±aj±Q are obtained by sampling (3 2) at x % ± aj ± a , 1 e , = 

t n ) Similarly to the one-dimensional case, the evolution points are m smooth 
regions and therefore the integral on the RHS of (3 3) can be replaced with a sufficiently 
accurate quadrature such as the Simpson rule (2 4), which leads to a scheme that is 
fourth-order accurate m time The derivatives at time t n , <p[± aj±Q are obtained by 
sampling the derivative of the reconstruction (3 2), i e , <p' t ± aj ± a = p'(x,±aj± Q , t n ) The 
other intermediate values of the derivative in time that are required m the quadrature can 
be predicted using a Taylor expansion or with a Runge-Kutta method m an analogous 
way to the one-dimensional case 

Remarks 

1 We present two different algorithms for constructing <p t ± a>J ± a : two-dimensional 
interpolants defined on two-dimensional stencils and a dimension-by-dimension 
approach We present both algorithms for the third-order scheme and extend the 
simpler dimension-by-dimension approach to fifth-order Our numerical 
simulations in §4 indicate that both reconstructions of <p,±a,j± Q are of a 
comparable quality In both approaches, the reconstruction of the derivatives 
V<p t ± Q j±a is done dimension-by-dimension 

2 We reproject v^+Jj+a and back onto the integer grid-pomts, obtaining 

cp"^ 1 We present several ways to carry out this reprojection a genuinely 
two-dimensional approach, a dimension- by-dimension strategy and a reprojection 
along the diagonal line through x t - aJ - a and x l+aj+a 
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Figure 3.1: The location of the evolution points Xi± a j± a and the domain of definition 
of the interpolants (pi± a ,j± a in two dimensions. 

3.2 Two-Dimensional Third-Order Schemes 

In order to obtain a third-order scheme, we need a fourth-order reconstruction of the 
point- values at the evolution points r, ±aj ± tt . 

3.2.1 A two-dimensional reconstruction of (pi± a j± a 

In this section we present a two-dimensional fourth-order reconstruction of the point- 
values <Pi± a ,j±a- In principle, a two-dimensional cubic interpolant would provide a recon- 
struction with the desired accuracy. Such an interpolant is based on a ten-point stencil. 
As usual, solving such a direct interpolation problem is unsatisfactory as spurious os- 
cillations might develop as a result of the lack of smoothness in the solution. Instead, 
we generate a two-dimensional fourth-order reconstruction as a convex combination of 
four quadratic interpolants, each which is based on a six-point stencil. We choose com- 
pact quadratic interpolants such that the union of all the six-point stencils is a compact 
ten-point stencil. Similarly to any WENO-type reconstruction, when singularities are 
present the six-point stencils containing the singularities are suppressed. In any case, 
we implicitly assume that the solution is sufficiently resolved such that the singulari- 
ties in the solution are isolated in the sense that they do not occur along neighboring 
parallel cell edges. Singularities will in general occur along adjacent cell edges. There 
is a lot of flexibility in choosing the ten-point stencil as well as the different six-point 
stencils. Here, for the evolution point x t+QJ+ „ we choose the ten-point stencil shown in 
Figure 3.2. We choose to use the four six-point stencils that are shown in Figure 3.3. 
Obviously, the union of these stencils is the ten-point stencil in Figure 3.2. Furthermore, 
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j+2 

j+i 

j 

j-i 


Figure 3.2: The ten-point stencil for the two-dimensional reconstruction of ^+ aj -+ a . 
The open circle shows the location of the evolution point at Xi +a j +a . 

they all enclose the cell containing the evolution point and they all cross different edges 
of the enclosing cell. A singularity along an edge will suppress two of these stencils, 
while a singularity in a corner will suppress three of these stencils. 

Remarks . 

1. The stencils for the evolution point at Xi_ Q j_ Q are obtained by a rotation of 180 
degrees of the stencils in Figures 3. 2-3. 3. 

2. We could use less than four stencils and still generate a scheme that will have the 
desired order of accuracy. 



i-i i i+1 i+2 


Given the four six-point stencils in Figure 3.3, a straightforward computation shows that 
third-order approximations for smooth <p at the evolution points Xi± a j ± a , <p$± a j± a = 
(p Vj±a) + O ( Ax 3 , Ay 3 ), Wk €{1,2, 3, 4}, is obtained with 


( Pi±a t j±a ~ a l Vhj "I" 
( Pi±aJ±a = GsViJ + 
( Pi±a,j±a = <PiJ + 
( Pi±aj±a = a 5 <PiJ + 

where 

Ci\ = 1 — 3a ■+■ 2 a: 2 , 

1 1 2 

a, = -5“+2“. 

a 7 = 1 — 2a 4- a 2 , 


a 2 ( Pm,j + a,2<Pij±i + a3<pi±ij±i 
a 6 ( Pi±l,j + &2 l Pi,j±l + a 3<A:fcl,j±l 
a 2<Pi±l,j + 0-2^Pij±l + 08^t±l, j±l 

®2^i±l,j 4"" 1 4” 


+ 0>4<Pi±2j + 0’4 ( Pi,j±2, (3-4) 

+ a 4 i Pi,j±2 + a4<Pi^l,j, 

+ a 4fi±l,j^l + a 4fi^l,j±l, 

+ a 4 i Pi±2,j + 1) 


a2 = 2a — 2a 2 , 03 = a 2 , 

, 3 1 2 1 1 2 

o 5 = 1 - -a + -a , a 6 = -a--a 

as = —a + 2a 2 . 


(3.5) 
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Figure 3.3: The four six-point stencils that cover the ten-point stencil for the two- 
dimensional reconstruction. 


The linear combination 

4 

^ = V (%i±Gn 2/j±a) “1“ O (Ax , A y ) , 

k = 1 

is fourth-order accurate provided that the constants c* are taken as 

ci = ^ (5a - 1) , c 2 = c 4 = | (-2a + 1) , c 3 = a. (3.6) 

A two-dimensional CWENO reconstruction is a straightforward generalization of the 
one-dimensional case (compare with (2. 7), (2. 8)), 

4 

/c=l 

V. 

Here 

k 

k ^i±a,j±Q fc Ck 

W i±aJ±cc ~ ^4 / 5 a t±QrJ±a “ /_ , Qk \P 5 

2^/=l ^:±a,j±a ' ^tdbaj'iaJ 

with the constants Cfc given by (3.6). As usual, the smoothness measure for every stencil 
is taken as a normalized sum of the discrete L 2 -norms of the derivatives. If we define 
the forward and backward differences A +<Pij = <Pi+ij — <Pij, A 
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A yfij = <Pij + 1 — <fij, Ayifij = ifiij — <Pij- 1, then the smoothness measures for the 
evolution point x l+aj+a are given by 

$i +a,j+a = (Ax¥>ij) 2 + (A+<Pi+l,j) +(A+tpi t j+i) + (Ay Ifii'j) + (A+yJjj+i) 

+ (Ay <Pi+i,j ) 2 + ~ [(A+AJw+u) 2 + (Ay A~ y>i,j + x) 2 ] , 

Sf +a , j+a = (A+ipi ,j ) 2 + (A tVi-uf + (A^i,i+i) 2 + (A + (A+wj+i) 2 

+ (Ay V’i+i.j ) 2 + ^ [(At A-^ j ) 2 + (A+A-^- +1 ) 2 ] , 

Sf+aJ+a = (Ajv’ij) 2 + (Ajv’t.J+l) 2 + (A+yJ.-l,j+l) 2 + (Aj^t.j) 2 + (A+lfii+ij) 

+ (Ay Ifii + lJ-l) 2 + ^2 [(A^AJ lfiiJ + \) 2 + (Ay A“l/? i+ lj) 2 ] , 

S i+a,j + a = (Ajv>ij) 2 + (A+^i + lj) 2 + (A t<PiJ+l) 2 + (Ay V’i.j) 2 + (Ay 
+ (Ay¥>i+ij) 2 + ^2 [(A* AJv’i+ij) 2 + (Ay Ay V’i.j) 2 ] • 

The smoothness measures for the evolution point x,_ Qi j_ a are 

Si-a,j-a — (Ajv>i— 2,j) 2 + (Ajv>i_i,j) 2 + (Aj^i_i t3 _i) + (A+V’t,j-2) + ( Ay V’i.j- 1 ) 

+ (A+^-u-0 2 + ^ [(A+A^-u) 2 + (Ay A~ v’t.j-i) 2 ] , 
s i- a ,j-a = (A+v>i ,>) 2 + (Aj^i-1,3') 2 + (A+VJi-ij-i) 2 + (A+v>i-u) 2 + (Ay v’i-i.i-i) 2 

+ (Ay^i,i-2) 2 + ^2 [(AjA“v>»j) 2 + (Ay Ay v’m-i) 2 ] , 

S i-a,j-a ~ (Aj^i-lj) 2 + (A+Vij-l) 2 + (A+<Pi-l,j-l) 2 + (Ay (fij-i) + (Ay^i-lj) 

+ (Ay Vi-ij-i) 2 + ^2 [(AjAjv»ij-i) 2 + (Ay Ay i,j) 2 ] > 

s i-a,j-c = (A+v»<— 2j) 2 + (A+Vi-1,,-) 2 + (Aj^i_ij_i) 2 + (A+^j) 2 + (A+^j_i) 2 

+ (A+^-u-x ) 2 + ^ [(A+AJ^-m ) 2 + (Ay A“ <Pi,j) 2 ] • 


3.2.2 A dimension-by-dimension reconstruction of <Pi± a ,j±a 

A different way to obtain high-order approximations for the values of <p l ± a j± a is by 
carrying out a sequence of one-dimensional reconstructions from §2.2. This dimension- 
by-dimension approach for the reconstruction step is similar in spirit to that of [17], but 
here, in order to generate a Godunov-type scheme (unlike [17]), we are forced to use 
evolution points that are not positioned in the same locations as the data x tJ . An ap- 
propriately chosen sequence of one-dimensional reconstructions addresses this problem. 

We use the subscript V to denote the full range of an array, such that <p,j and 
<Pi t * denote the one-dimensional arrays <p*j = (<fiij, . . . , <Pnj) and = (<A,i, • • • , <A,aO- 
With the notation for the one-dimensional third-order reconstruction, (2.10), we can 
express the dimension-by- dimension reconstruction at x i+a j +a as 

1. For each i,j: <fii+ a ,j = reconstruct_cp_lD_3 (i, a, <p t j) 

2. For each i,j: '~p l+a ^ +a — reconstruct_c^_lD_3 (i, o, <£»+<*,»)• 
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Figure 3.4. The dimension-by-dimension reconstruction process m two dimensions Left: 
the first step where the intermediate interpolants p %+a ,] at x l+aj (open squares) are 
computed using the data <p tiJ (black dots) Right the second step, where <p t + a ,j is 
interpolated in the j direction, giving <^ t+Q)J+a at x l+a ^ J+a (open circle) 

Here, we first interpolate along the first coordinate axis and reconstruct p at x l+aj 
The data at x t+oc<J is then interpolated along the second coordinate axis to the location 
x i+aj+a to give ipi + aj +Q (see Figure 3.4) Obviously, the order m which the steps are 
performed is not important In a similar way, a dimension-by-dimension reconstruction 
at rr t _ Q)J _ a is given by 

1 For each i , j ip t - a j = reconstruct.^ _1D_3 (i, — a , P*,j) 

2 For each i, j p t - aj - a — reconstruct.^.l D .3 (i, —a, <p x - a ,*) 

3.2.3 The reprojection step 

After evolving the solution to the next time step at the evolution points we 

would like to reproject p'f+a J+a back onto the integer grid points x h] to end up with 
There are several different ways to perform this task out of which we choose to present 
the following a two-dimensional reprojection using the two-dimensional reconstruction 
of §3 2 I or the dimension-by-dimension reconstruction of §3 2 2, and a one-dimensional 
projection along the diagonal 

1 A 2D reprojection. The evolution points at aj,± aj ± a have the same geometrical 
relationship to x X] as x tJ has to Hence, m order to reconstruct ip”* 1 

from <Pi±aj±cn we can directly utilize the projections from §3 2 1 or §3 2 2, taking 
<p t ±a,j±a as the input data, and reversing the sign of the parameter from ±a to 
The final value is then taken as the average of the projections of <p t+QJ+a 
and Hence, if we denote either the two-dimensional or the dimension- 

by-dimension reconstruction described in §3.2 1 or §3 2 2 as 


c P%±aj±a = reconstruct_<p^D_3 («, j, ±a, <p ) , 


(3 7) 


where tp is now the two-dimensional array then the reprojection step is 
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(a) For each i,j = reconstruct j/p_ 2D_3 (z, —a, ^+ aj+a ) 

(b) For each i,j = reconstruct _3 (*, a, <p t - a j-a ) 

(c) For each i,j tf* 1 = | (y>+ + <?“ ) 

2 A diagonal reprojection. In this case we use one-dimensional data along the di- 
agonal, ■{</?»— 1 + 0 ,^— i+aj a,j— ai Pi+aj+ai ^Pi+i— a,j+i— a}) construct a third-order 

WENO approximation of ip™* 1 (see Figure 3 5) 

Define 


or a — 1 

2a — ^-i+ Q>J -i + a + 2{2a-l) lp '- aJ - e 


1 _ cx ■" _ - 

H ~ ( Pi+aj+a ~ (^ij) + O (Ax , A?/ ) , 


1 — a 


Of — 1 


2 2 (2a-l) Vt+aj+a 

q : 2 

+ g— j^+l-aj+l-a = Kj) + O (Ax 3 , Ay 3 ) 


Since (</? tJ + <^)/2 = (/? (x tJ ) + O (Ax 4 , Ay 4 ), we can obtain <p "+ 1 as 




(3 9) 


where as usual /(a^+a ) > and a to = (2 (c + S^) P ) 1 The smoothness 

measures are again taken as the sum of the discrete L 2 norm of the derivatives, 
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which in this case is more complicated due to the uneven spacing of the data 
1 


^ “ Ax 


Ax 3 


) 2 / \ 2 
. / tPi+aj+a (Pi—aj—i 

v 2 a 

4 f ct,j— a Cft-l+qj-l+q ^i+qj+tt ‘ft— aj— a 

1 -2a 2a 


5+ = -f- 
J Ax 


( ft+a.j+a ft-a,j-a 

V 2a 


ft+l-aj-H-a ft+aj-f a A 


+ 


1 -2a 


J 


4 f ft+aj+a ft— a,j- a ft+1— aj+l— a ft+a,j+a 

Ax 3 \ 2a 1 — 2a 


Remark Our numerical simulations in §4 3 indicate that there is little difference between 
the quality of the two-dimensional reconstruction and the dimension-by-dimension re- 
construction of §3 2 1 and §3 2 2 We will use this fact when extending our methods to 
fifth-order and higher dimensions We note that the diagonal reprojection significantly 
reduces the CFL number (see §4 4) 


3.3 A Two-Dimensional Fifth-Order Scheme 

Using the dimension-by-dimension approach, it is easy to extend the above scheme to 
fifth-order simply replace the one-dimensional third-order interpolations by the fifth- 
order interpolation in §3 2 2 Using the one-dimensional notation, (2 19), we obtain a 
fifth-order reconstruction at x l+a>J+a as 

1. For each z,j ip t+a>J = reconstruct _</>_lD_5 (z, a, 

2. For each t, j (p i+a>J+a = reconstruct j^_1D_5 ( j , a, <Pi+ a ,*) 

Similarly, at we have 

1. For each z, j = reconstruct _<^_1D_5 (z, —a, 

" ' 2. For each z, j ~<p t - a ,]-a ~= recoristruct‘_^_l‘D_5 ipr~~a]*) 

We denote this reconstruction as 

ft±a,j±a = reconstruct.<p_2D_5 (z,j, ±a, <p) (3 10) 

For the derivatives we have 

1 For each z, j (p' l±Qt} = reconstruct V-1D-5 (z, ±a, 

2. For each z, j ip' l±aj±ac = reconstruct _y?'_lD_5 (z, ±a, ip x ± a ,*) 

which we denote as 

^± tt ,j±Q = reconstruct V-2D.5 (z, j, ±a, <p) (3.11) 

Reprojection onto the original grid points x tJ is performed using the two-dimensional 
dimension-by-dimension reprojection option described m §3 2 3 
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Remarks 

1 Due to the reduced stability resulting from the use of diagonal reprojection, 
which is demonstrated in §4.4, we do not develop a fifth-order analog to the 
third-order diagonal reprojection 

2 It is straightforward to develop a fifth-order two-dimensional method involving 
two-dimensional stencils, extending §3 2 1 Such a method would involve four 
interpolants defined on ten-point stencils that cover a twenty one point stencil. 

We summarize the two-dimensional fifth-order algorithm in the following, where 
RK At) is now the fourth-order RK method which integrates (2 1) As in 

Algorithm 2.1, each internal step of the RK method will require additional reconstruc- 
tions of (p[± a from that step’s (p t ± a 

Algorithm 3.1 Let a = 1/ (2 + \/2)* Assume that {<p" } are given 

1 Reconstruct 


<Pi±aj±a = reconstruct^ -2D .5 (i, j, ±a, <p) 


= reconstructs^ -2D-5(t,j,±a, cp) 


2 Integrate 


( Pi ± a i j±ot — RK (^±a,j±a> a,j±a> At/2) 

V’.iajia = reconstruct.^'.. 2D. 5 (t, ±a, v?±tf±<*) 

Vwfcajia = RK (S±a,j±a> L P'i±aj±a> At) 

^t±a]±a = reconstruct-^’ .ZD .5 (*, ±c*, <P±t|±a) 

= ¥>3^*. + Y [h m (vd*.) + U 


3. Reproject 


= reconstruct.<p.2D.5 (i,j, <^£*± 0 ) 
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3.4 Multi-Dimensional Extensions 

The extension of the direction-by-direction approach to more than two space dimensions 
is straightforward. For example, using the notation of §3 3, a three-dimensional fifth- 
order reconstruction is 

1 For each z, j, k. <p t+aiJifc = reconstruct .</? _1D_5 (z, a, 

2 For each i, j, k <^, +aj+a>fe = reconstruct.^_lD_5 (z, a, <p,+ a ,*,k)- 

3 For each i, j, k (p t + a ,j+a,k+a = reconstructj^.lD.5 (z, a, ^ +aj+a ,*) 

The reconstruction at x x _ a ] _ a ^^ a is handled similarly, and the same for the reconstruc- 
tion of <Pt+aj+a,fc+Q I n three dimensions a = 1/ (3 4- \/3) 

A d-dimensional reconstruction based on d-dimensional stencils quickly becomes very 
large It is readily apparent that the dimension- by- dimension approach will scale to high 
dimensions better than d-dimensional interpolants 


4 Numerical Simulations 

In this section we present simulations that test the schemes we developed in this paper 
In §4 1 we demonstrate the third- and fifth-order method in one dimension §4 2 focuses 
on the fifth-order method m two and three space dimensions In §4 3 we compare 
the two-dimensional third-order method based on two-dimensional stencils with the 
direction-by-direction approach In §4 4 we examine, m detail, stability issues in two 
dimensions, including comparisons with [17] Some of these examples are standard test 
cases that can be found, e g , in [22, 31, 35] 

We do not follow the practice m [17] of masking singular regions from our error 
measurements. 

4.1 One-Dimensional Examples 

A- convex- Hamiltonian 

We start by testing the performance of our schemes on a convex Hamiltonian. We 
approximate solutions of the one-dimensional equation 

4>t + ^ (4>x + l) 2 = 0, (4 1) 

subject to the initial data <p(x, 0) = —cos(rx) with periodic boundary conditions on 
[0,2] The change of variables, u(x,t) = cj> x (x,t) + 1, transforms the equation into 
the Burgers’ equation, u t + \ (u 2 ) x = 0, which can be easily solved via the method of 
characteristics [35] As is well known, Burgers’ equation generally develops discontinuous 
solutions even with smooth initial data, and hence we expect the solutions of (4.1) to 
have discontinuous derivatives In our case, the solution develops a singularity at time 

t = 7T -2 
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The results of our simulations are shown in Figure 4 1 The order of accuracy of 
these methods is determined from the relative L 1 error (see [30]), defined as the ZJ-norm 
of the error divided by the Z^-norm of the exact solution These results along with the 
relative L°°-norm before the singularity, at T = 0.8/n 2 , are given in Table 4 1, and after 
the singularity at T = 1 5 / 7 r 2 m Table 4 2 


N 

Third-order method 

relative L 1 -error 

L 1 -order 

relative L°° -error 

L°°-order 

100 

9.41 xl0“ 5 

- 

1 77xl0~ 5 



200 

1.13xl0 -5 

3 06 

1 33xl0 -e 

3 73 

400 

1.39xl0" 6 

3 02 

9 35x10"* ” 

3 83 

800 

1.74xl0' 7 

3 00 

5 94xl0" y 

3 00 

N 

Fifth-order method ”] 

relative L l -error 

L 1 -order 

relative L°° -error 

L°°-order 

100 

1 41xl0 -5 

- 

2 61 xl(r e 

- 

200 

4 21 xlO -7 

5 07 

4 03x10”* 

6 02 

400 

3 31xl0 -8 

5 00 

6 53xl0 -lu 

5 95 

800 

4 03xl0 _lu 

5 03 

1 OOxlO" 11 

6 03 


Table 4 1: Relative L 1 -errors for the one-dimensional convex HJ problem (4 1) before 
the singularity formation T = 0 8/x 2 . 


TV 

Third- order method 

relative L l -error 

L 1 - order 

relative L°° -error 

L°°-order 

100 

9 10xl0 -4 

- 

2 77xl0 -4 

- 

200 

2 16xl0 -4 

2 07 

7 63xl0 -5 

1 86 

400 

...6 84xl_0r 5 

166 

2 68xl0 _5 ___ 

__1 51 

800 

2 75xl0 -5 

1.31 

2.08xl0- 5 

0.37 

N 

Fifth-order method 

relative L l -error 

L L -order 

relative L°° -error 

L°°-order 

100 

7 85xl0" 4 

- 

5 78xl0 -4 

- 

200 

1 61xl0 -4 

2.29 

8 29x10"* 

2 29 

400 

6 71xl0 -5 

1 26 

5 09xl0- s 

1 26 

800 

3 44xl0 -5 

0 96 

3 44x10"^ 

0 96 


Table 4 2 Relative L 1 -errors for the one-dimensional convex HJ problem (4 1) after the 
singularity formation T — 1 5 /vr 2 
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Figure 4 1 One-dimensional convex Hamiltonian (4.1) Left: the solution before the 
singularity formation, T = 0 8/n 2 Right the solution after the singularity formation, 
T — 1 5 / 7 T 2 N = 40 Shown are the third- and fifth-order approximations, and the 
exact solution. 

A non-convex Hamiltonian 

In this example we deal with non-convex Hamilton-Jacobi equations In one dimension 
we solve 

cf>t - cos (4> x + 1) = 0, (4 2) 

subject to the initial data <j> ( x , 0 ) = — cos (nx) with periodic boundary conditions on 
[ 0 , 2 ] In this case (4 2) has a smooth solution for t < 1 049/7T 2 , after which a sin- 
gularity forms A second singularity forms at t « 1 29/7 r 2 The results are shown in 
Figures 4 2 The convergence results before and after the singularity formation are given 
in Tables 4 3-4 4. 


A linear advection equation 


In this- example -(-[1 7] with- a misprint, corrected in [40]) we solve the one-dimensional 
linear advection equation, 1 e , H (<p x ) = <p x We assume periodic boundary conditions 
on [—1,1], and take the initial data as (x, 0) = g (x — 0 5) on [— 1 , 1 ], where 


$(*) = “ 


'y/Z 9 2tt' 

~{- — -f- 

2 2 3 


(2 + 1 ) + h(x), 


h(x) = 


2 cos (y x2 ) ~ A 
3/2 + 3 cos (27 tx) , 

15/2 — 3 cos (2ttx) , 

(28 + 47 r + cos (3ttx)) /3 + 6nx ( x 


-1 < x < 

— 5 < x < 0 , 
0 < x < |, 

1 ) , 1 < x < 1 


(4 3) 


The results of the fifth-order method are shown in Figure 4 3, where it is compared with 
the fifth-order method of [17] The reduced dissipation effects of our method are visible 
in the reduced round-off of the corners 
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N 

Third-order method 

relative L 1 -error 

L l -order 

relative L°° -error 

L°°-order 

100 

6 47xl0 -5 

~ 

9 05xl0 -6 

- 

200 

7 78x10“® 

3 06 

1 llxl0 -ti 

3 03 

400 

8 77xl0" 7 

3 15 

9 27xl0 -8 

3 58 

800 

9 87xl0 -8 

3 15 

6 12xl0 -9 

3.92 

N 

Fifth-order method 

relative L 1 - error 

L l -order 

relative L°° -error 

L°°-order 

100 

1 29xl0“ 5 

- 

4 97xl0“ 6 

- 

200 

6 52xl0“ 7 

4.31 

2 38xl0 -7 

4 38 

400 

2 10xl0 -8 

4.95 

6 13xl0" 9 

5 28 

800 

5 96xlO“ 10 

5 14 

1 03x 10~ 10 

5 90 


Table 4 3: Relative ZJ-errors for the one-dimensional non-convex HJ problem (4 2) 
before the singularity formation T = 0 8/7T 2 


N 

Third-order method 

relative L 1 -error 

L l -order 

relative L 00 -error 

L°°-order 

100 

2 81xl0 -4 

- 

9 64xl0 -5 

- 

200 

1 32xl0 -4 

1 08 

5 05xl0 -5 

0 93 

400 

2 31 xlO -5 

2 52 

6 OOxlO -6 

3 07 

800 

8 43xl0 -6 

1 46 

3 30xl0~ 6 

0 86 

N 

Fifth-order method 

relative L 1 -error 

L 1 -order 

relative L°° -error 

L°°-order 

100 

1 57xl0" 4 

_ 

1.12xl0~ 4 

- 

200 

8 34xl0 -5 

0.91 

6 60xl0~ 5 

0 77 

400 

1 22xl0 -5 

2 78 

8 64x10-® 

2 93 

800 

6 67xl0 -5 

0.87 

5 23x10-® ' 

072 


Table 4 4 Relative ZJ-errors for the one-dimensional non-convex HJ problem (4 2) after 
the singularity formation T = 1.5/ 7r 2 
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Figure 4 2 One-dimensional non-convex Hamiltonian (4 2) Left The solution before 
the singularity formation, T = 0 8 / 7 r 2 . Right. The solution after the singularity forma- 
tion, T = 1 5/tt 2 . TV = 40 Shown are the third- and fifth-order approximations, and 
the exact solution 


t = 2 t = 8 



Figure 4 3 One-dimensional linear advection, (4.3) T = 2, 8, 16, 32. TV = 100 Crosses 
our fifth-order method Circles the fifth-order method of [17] with a local Lax- Friedrichs 
flux Solid line the exact solution 
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4.2 Two-Dimensional Examples 


A convex Hamiltonian 

In two dimensions we solve a problem similar to (4 1) 


+ fiy + l) 2 — - 0 , 


(4 4) 


which can be reduced to a one-dimensional problem via the coordinate transforma- 


tion 


1 1 
1 -1 


The results of the second-order calculations for 


the initial data <£(x, y, 0) = — cos (7r(x + y)/2) = — cos ( 7 r£) are shown m Figure 4 4 
The convergence rates for the two-dimensional fifth-order scheme before and after the 
singularity are shown in Table 4.5 


Before singularity T = 0 8/7T 2 

N 

relative L l -error 

L l -order 

relative L°° -error 

L°°-order 

50 

1.19xl0“ 4 

- 


Hi 

■ESI 


4 13 


5 56 

JEM 

I 

5 30 

1 

7 20 

After singularity T — 1.5/7T 2 

N 

relative L l -error 

L l -order 

relative L°° -error 


50 

1.32xl0~ 3 

- 

2.07x10“* 

- 

100 

3.89xl0 -4 

1 76 

3 60x10"* 

2 52 

200 

4.86xl0 -5 

3 00 

1 69xl0 -7 

4 41 


Table 4 5 Relative L 1 - and Z/°°-errors for the two-dimensional convex HJ problem (4 4) 
before and after singularity formation, computed via the fifth-order method 


A non-convex Hamiltonian 

The two-dimensional non-convex problem, which is analogous to the one-dimensional 
problem7(4" 2)7 is 

<f>t - cos (<j i) x + <f) y + 1) = 0. (4 5) 

Here we assume initial data that is given by <p (x, y, 0) = — cos ( 7 r(x + y ) /2), and periodic 
boundary conditions The results are shown m Figure 4 5 The convergence results for 
the two-dimensional fifth-order scheme before and after the singularity formation are 
given m Table 4 6 

A fully two-dimensional example 

The above two-dimensional examples are actually one-dimensional along the diagonal 
To check the performance of our methods on fully two-dimensional problems we solve 

<j>t + 0X02/ = o, (4 6) 
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2 ? 


2 2 


Figure 4.4 Two-dimensional convex Hamiltonian, (4 4) Left the solution before the 
singularity formation, T = 0 8/tt 2 . Right: the solution after the singularity formation, 
T = 1 5/-7T 2 N = 40 x 40. The solution is computed with the fifth-order method 


Before singularity T — 0 8/7r 2 

N 

relative L l -error 

L l -order 

relative L°° -error 

L°° -order 

50 

1 llxlO -4 

- 

1 26xl0~ 6 

- 

100 

6 91xl0 -6 

4.00 

2 42xl0 -s 

5 70 

200 

3 85xl0~ 7 

* *4.17 ~ 

“ 6 27xl0 -lu 

' 5 27' 

After singularity T — 1 5/n 2 

N 

relative L 1 -error 

L l -order 

relative L°° -error 

L x '-order 

50 

1 47xl0 -3 

- 

8 58xl0 _s 

~ 

100 

1 93xl0' 4 

2.93 

9 27xl0 -7 

3 21 

200 

8 87xl0 -5 

1.12 

3 09xl0 -7 

1 58 


Table 4 6 Relative L l - and L°° -errors for the two-dimensional non-convex HJ problem 
(4 5) before and after the singularity formation, computed with the fifth-order method 
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Figure 4 5 Two-dimensional non-convex Hamiltonian, (4 4) Left, the solution before 
the singularity formation, T = 0 8/7T 2 Right the solution after the singularity for- 
mation, T = 1.5/ 7r 2 N = 40 x 40 The solution is computed with the fifth-order 
method 

on [ — 7r, 7r] x [— 7r,7r], subject to the initial data <p(x,y,0) = sin(x) -1- cos (y) with pe- 
riodic boundary conditions The exact solution for this problem is given implicitly by 
4> (x, y,t) = — cos (i?) sin (r) + sm (q) + cos ( r ) where x = q-t sin (r) and y = r + t cos (q). 
This solution is smooth for t < 1, continuous for all t and has discontinuous derivatives 
for t > 1 The results of our simulations at times T = 0 8, 1 5, are shown in Figure 4 6 
The convergence results for the first- and second-order two-dimensional schemes be- 
fore the singularity formation are given in Table 4 7 and confirm the expected order of 
accuracy of our methods 


Before singularity T — 0 8 

N 

relative L 1 -error 

L l -order 

relative L°°-error 

L°° -order 

- 50 

6 10x10-°- 

----- 

8 15xl0 -s 

~ 

100 

2 10xl0" 7 

4 86 

7 35xl0~ lu 

6 79 

200 

7 53X10- 9 

4 80 

5 59x10” 12 

7 04 


Table 4 7 Relative iJ-errors for the two-dimensional HJ problem (4 6) before singularity 
formation T = 0 8 The solution is computed with the fifth-order method 


An eikonal equation in geometric optics 

We consider a two-dimensional non-convex problem that arises in geometric optics [20] 

(4 7) 


f 4>t + ~ 0’ 

\ ( f > (x, y, 0) = \ (cos (2ttx) - 1) (cos (27 ry) - 1) - 1 
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Figure 4 6 Fully two-dimensional Hamiltonian, (4 6) Left the solution before the 
singularity formation, T = 0 8 Right the solution after the singularity formation, 

T = 15 N = 50 x 50 The solution is computed with the fifth-order method 

The results of our fifth-order method at time T = 0 6 are shown m Figure 4 7, where 
we see the sharp corners that develop in this problem 

An optimal control problem 

We solve an optimal control problem related to cost determination [35] Here the Hamil- 
tonian is of the form H(x, y, V(j > ) 

/ 4>t ~ sin (y) (j> x + sin (: x ) <j> y + \cf> y \ - ± sin 2 (y) - 1 + cos ( x ) = 0, , , 

- I M*,v, 0) = 0 . ..... . _ .. , (48) 

The result of our fifth-order scheme is presented m Figure 4 8 and is in qualitative 
agreement with [31] 

4.3 A Comparison of Two-Dimensional Third-Order Interpolants 

In this section we use the examples (4.4), (4 5) and (4 6) to compare the third-order 
method of §3 2 1, based on interpolation via two-dimensional stencils, with that of §3 2 2, 
where we used a direction-by-direction approach The results are shown in Table 4 8 
The dimension-by-dimension method produces errors that are approximately twice as 
large compared with the genuinely two-dimensional reconstrcution However, the con- 
vergence rate is qualitatively the same in both methods These results motivated us to 
base our fifth-order scheme on the much simpler dimension-by-dimension reconstruction 


High Order Schemes for HJ Equations 


31 




Figure 4 7 Two-dimensional eikonal equation, (4 7) N = 40 x 40. Left the initial 
data Right: the fifth-order approximation at T = 0 6 



4 -4 


Figure 4 8 Two-dimensional optimal control problem, (4 8) An approximation with 
the fifth-order method is shown at T = 1 N — 40 x 40 
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N 

2D Stencils 

Direction-by- Direction 

relative L l -error 

L l -order 

relative L 1 -error 

L l -order 

Convex Hamiltonian atT = 0 8/n' 2 

50 

4 70xl0 -4 

- 

6 13xl0 -4 

- 

100 

7 54x10“* 

2.64 

9 43xl0 -5 

2 70 

200 

8 07xl0 -6 

3.23 

1 02xl0 -5 

3 21 

Convex Harmltoman atT = 1 5/7r 2 

50 

1 23xl0 -3 

- 

2 61xl0" 3 

- 

100 

4 56xl0 -4 

1 44 

8 19xl0 -4 

1.67 

200 

3 70xl0~ 5 

3 62 

1 22xl0" 4 

2.74 

Non- Convex Hamiltonian atT — 0 8/n 2 

50 

2 27xl0 -4 

- 

3 92xl0“ 4 

- 

100 

3 75xl0 _& 

2 60 

6 97xl0 -5 

2 49 

200 

3 99xl0 -e | 

3 23 

7 22xl0 -e 

3 27 

Non- Convex Hamiltonian atT = \5/n 2 

50 

1 23xl0 -3 

~ 

1 94xl0 -3 

- 

100 

2 50xl0 -4 1 

2 30 

4 16xl0~ 4 

2.22 

200 

7 63xl0 -5 

1 71 

1 20xl0 -4 

1 79 

Fully 2D Example atT = 0 8 

50 

2 OlxlO -4 

- 

1 48xl0“ 4 

- 

100 

2 42x 10 -5 

3 05 

1 65xl0" 5 

3.16 

.200 

~2 95xl0- 6 ' 

3 04 . 

1J5x10"1_ 

3 08 


Table 4.8 Comparison of the third-order method of §3 2 1, with an interpolation via 
two-dimensional stencils, and that of §3.2.2, with the direction-by-direction approach 
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4.4 A Stability Study 

In this section we present a couple of stability studies we obtained in our simulations 
We start by checking the stability properties of the third-order scheme with different 
reprojection steps The reconstruction step is done m all cases using the direction-by- 
direction interpolant We compare the dimension-by-dimension reprojection and the 
diagonal reprojection (of §3 2 3) In Figure 4 9 we plot the L 1 error as a function of the 
CFL number. The test problem is (4 6) with the fully two-dimensional Hamiltonian 
The solution is computed at T = 0 8 We see that the use of a diagonal reprojection 
significantly reduces the maximum allowed CFL number 


Third order fully 2D H, T=0 8 



Figure 4 9: Stability of the two-dimensional third-order method with a dimension-by- 
dimension ( crosses ) vs a diagonal reprojection ( diamonds ) Fully two-dimensional 

Hamiltonian (4.6 )- -T =08 (before singularity) N =400 x 100- - 

We now turn to checking the stability properties of the two-dimensional fifth-order 
method of §3 3 by computing the L 1 errors for various examples while varying the CFL 
number. In Figure 4 10 we compare the results obtained with our fifth-order scheme 
with the fifth-order method of [17], for which we used a local Lax-Friedrichs flux The 
numerical tests indicate that larger CFL numbers can be used with our method 

4.5 Three-Dimensional Examples 

We proceed with a three-dimensional generalization of the convex Hamiltonian (4 4), 
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x1 q- 4 non-convex H, T=0 8/jc 2 xIO" 3 non-conve* H, T=1 5/n 2 




Figure 4 10: Stability of the two-dimensional fifth-order method N = 100 x 100. 
Crosses our fifth-order method Circles the fifth-order method of [17] with a local 
Lax-Friedrichs flux Upper left linear advection ( H (V</?) = Vip) with initial condition 
<f> {x. y, 0) = — cos (tt(x + y)/ 2) Upper right fully 2D Hamiltonian (4 6) Middle row 
convex Hamiltonian (4 4), before the singularity (right) and after the singularity (left) 
Bottom row non-convex Hamiltonian (4 5), before the singularity (right) and after the 
singularity (left) 
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subject to the initial data <fi {%> V, z, 0) = — cos {n{x + y + z) /3) The convergence results 
for the three-dimensional fifth-order scheme before and after the singularity formation 
are given in Table 4 9. We also approximate the solution of the non-covex problem 

(fit ~ cos {(fi x + <fi y -f- (fi z + 1) = 0, (4 10) 

with the same initial data The convergence rates for the three-dimensional fifth-order 
schemes are given in Table 4 10 


Before singularity T = 0.5/7T 2 

N 

relative L 1 -error 

L l -order 

relative L°° -error 

L 00 -order 

25 

2 61xl0” 4 


1.07xl0" 7 

~ 

50 

6 40x10”“ 

5 35 

3 16xl0" lu 

8 41 

100 

1 50xl0” 7 

5 42 

9 18xl0 -13 

8 43 

After singularity T = 1 5/7 r 2 

N 

relative L l -error 

L l -order 

relative L°° -error 

L°° -order 

25 

6.95xl0” 3 

- 

1 80xl0” 5 

- 

50 

1 40xl0” 3 

2 31 

4.15x10”“ 

2 12 

100 

5 33xl0 -4 

1 39 

6.94xl0” 7 

2 58 


Table 4 9 Relative L 1 - and L°°-errors for the three-dimensional convex HJ problem 
(4 9) before and after the singularity formation, computed with the fifth-order method 


Before singularity T = 0 5/7 r 2 

N 

relative L 1 -error 

L l -order 

relative L°° -error 

L°° -order 

25 

7.28xl0 -4 

- 

3 70xl0” 7 

- 

50 

3.71X10" 5 

4 29 

4.06xl0" y 

6 51 

100 

1 05x10”“ 

5 14 

2 18x10”“ 

7 54 

After singularity T — 1 5/7T 2 

N 

relative L l -error 

L l -order 

relative L°° -error 

L 00 -order 

25 

6 74x10-*- 

- = = 

-3.27x10-“ 

~ 

50 

1.26x10-* 

2 42 

6.90xl0 -7 

2 25 

100 

4.21xl0” 4 

1 59 

6.84x10"“ 

3 33 ! 


Table 4 10. Relative L 1 - and L°°-errors for the three-dimensional non-convex HJ problem 
(4 10) before and after the singularity formation, computed with the fifth-order method 
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